A CONSERVATION SCIENCE PARTNERS REPORT
A summary of data products created for the US Forest Service, including links to the data and web application used to render deliverables, metadata regarding the accuracy and precision of specific results, and guidance for appropriate uses of the data.
The Four Forest Restoration Initiative (4FRI) is a landscape-scale project focused on restoring Ponderosa pine forests in Arizona. An important part of the successful implementation of this project is to assess the impacts of restoration treatments on forest structure. Northern Arizona University (NAU) and the United States Department of Agriculture (USDA) Forest Service collaborated to quantify and describe the amount, pattern, and distribution of canopy cover within and around Ponderosa pine forests of the South Kaibab and Coconino National Forests.
The project area encompasses 1,224,900 acres and intersects portions of both Coconino and Kaibab National Forests (Figure 1).
Figure 1: The project area / orthoimagery acquisition boundary in relation to the administrative boundaries of Coconino National Forest and the South Kaibab (Tusayan and Williams Ranger Districts). Selectable layers include image tile boundaries (n = 619) and the locations of the spatially balanced random survey cells used to develop training data for the canopy cover classification model (n = 158, ~300-acre cells). All told, 6,119 samples were used to train the model, 40% of which were collected completely at random (from cells in blue) and 50% of which were collected opportunistically (from cells in red). The remaining 10% of samples were collected completely opportunistically outside of the spatially balanced survey cells.
High quality aerial imagery was acquired for the project area between June 6 and June 23, 2014 by the USDA Farm Service Agency Aerial Photography Field Office (APFO). The acquisition platform was a light aircraft flying ~5,570 m above ground level to achieve a nominal resolution of 0.3 m. The 4-band – red, green, blue, and near infrared (hereafter, ‘NIR’) – imagery was collected using a Microsoft UltraCam Eagle sensor with a 100.5 mm focal length. The images were orthorectified, mosaicked, and radiometrically adjusted to meet contract requirements. The primary difference between this imagery and the better-known NAIP image archive is quality: i.e., how far off nadir acquisitions are allowed to be, the increased overlap between photos, a requirement for no cloud cover, and restrictions on time of day (to reduce shadows).
The APFO imagery was stored as an asset (i.e., a stack of images, referred to as an ImageCollection) in Google Earth Engine (hereafter, ‘Earth Engine’). Each ee.Image object within the collection represents a single Digital Orthophoto Quarter Quarter Quadrangle (DOQQQ).
We developed a 3-class (tree, non-tree, and shadow) supervized classification model. Specifically, we used a random forest classifier (Breiman 2001). The model-building process entailed several steps, each of which are described in detail below:
Spatially balanced random survey cells generated using the Reversed Randomized Quadrant-Recursive Raster algorithm (RRQRR; Theobald et al. 2007) were used to develop training data for the canopy cover classification model (n = 158, ~300-acre cells).1 All told, 6,119 samples were used to train the model, 40% of which were collected completely at random and 50% of which were collected opportunistically (Figure 1). The remaining 10% of samples were collected completely opportunistically outside of the spatially balanced survey cells.
A large suite of predictors were developed using the imagery as well as digital elevation data. For example, we computed the Normalized Difference Vegetation Index (NDVI) from the red and NIR bands in the imagery, and topographic layers (i.e., elevation, slope, aspect) using the USGS National Elevation Dataset (NED; Farr et al. 2007) digital elevation data. Additionally, we applied:
Figure 2: Examples of quantities derived using the imagery seen in true-color in the left-most panel. Predictors shown here (the singleband pseudocolor images in the right-most panel) include NDVI, DOG, entropy, and GLCM cluster shade.
To extract predictor variable information to the locations of all samples in the training data, we used reducers (i.e., the ee.Reducer class) in Earth Engine. Samples that fell in the overlapping area between DOQQQ tiles had two (or more) sets of covariate information. We did not allow these redundant copies enter the model-training step (described below). Instead, we selected only one set of covariate information for these samples, specifically the set corresponding to the DOQQQ tile whose centroid was nearest the ‘offending’ sample.
The training samples and associated covariate information were brought into R, where we used the caret package to find optimal parameters for the random forest algorithm. Specifically, we conducted a grid search of parameters, using repeated cross validation and accuracy as the performance metric, to select the optimal model.
The final, tuned model hyperparameters (the number of trees, variables per split, and minimum leaf population), were then used in the ee.Classifier.randomForest method in Earth Engine. We trained the model against a regionally-boosted sample for each image in the final prediction in order to better-calibrate the model against local conditions. The final canopy cover classification model can be viewed using the web application linked to below.
[TODO: insert screenshot of the web app!]
The performance of the classification model was evaluated against a ‘test’ partition, a set (n = 621) of samples selected at random from among the spatially-balanced training set and withheld from the model during tuning/training. Statistical measures of the performance of the model are provided in a confusion matrix (Table 1). Numbers along the diagonal (boldface font) represent correctly classified test samples. Overall accuracy (the sum of correctly classified samples divided by the total number of samples) was 96.5%.
Table 1: Confusion matrix for the classification model.| Predicted | |||
|---|---|---|---|
| Canopy | Shadow | Other | |
| Actual | |||
| Canopy | 195 | 1 | 11 |
| Shadow | 8 | 198 | 1 |
| Other | 1 | 0 | 206 |
Off-diagonal elements represent different types of errors. For example, there were 11 samples that were misclassified as ‘other’ (non-tree/non-shadow) when the test data show they were actually canopy. Additional measures of the performance of the classifier for each class are reported in Table 2. For example, sensitivity measures the proportion of the actual samples in a given class that were correctly identified as such, while the positive predictive value (or precision) is the proportion of predictions in a given class that were correct. For more information regarding performance measures in Table 2, see Wikipedia.
Table 2: Statistical measures of the performance of the model for each class. Class-wise statistics were computed using a ‘one against all’ approach.
| Sensitivity | Specificity | Positive predictive value | Negative predictive value | |
|---|---|---|---|---|
| Canopy | 0.942 | 0.978 | 0.956 | 0.971 |
| Shadow | 0.957 | 0.998 | 0.995 | 0.979 |
| Other | 0.995 | 0.971 | 0.945 | 0.998 |
Finally, we performed a ‘straight-face’ test (qualitative visual assessment) of the result. Though every iteration of the model we produced ultimately passed the statistical tests, we noticed that some regions within the project area were more problematic than others. For example, the craters in the NE quadrant of the study area were showing up with more area classified as tree canopy than expected, and a quick visual assessment confirmed that the model was likely confusing green grass in the understory as significant and contiguous enough to be canopy. This is what lead to the deployment of a larger opportunistic sample and the subsequent development of a geographically boosted training set.
Indicators of desired forest structural conditions (Science and Monitoring Working Group 2012) include patch size, density, and configuration. These landscape pattern indices (LPIs), among others, are described in Table 1. at multiple spatial scales.
The US Forest Service selected eight metrics (some of which were composites of several individual FRAGSTATS metrics) from Cushman et al. (2008).
We used the R package SDMTools to compute the majority of the LPIs.
The landscape/classification was divided up into sublandscapes. Landscape pattern indices (LPIs) were computed for each sublandscape and subsequently mosaiced into a new raster. LPIs were developed at multiple scales.
We developed the following mutually agreed upon spatial metrics for each area. TABLE. Accuracy and precision of results. Spatial scales. Very fine scale (1 acre) and a very very broad scale (watersheds, HUC5 or so) — 500 acres or so. Patches…. They (the people doing the treatements) allow for patches to be up to 4 acres. The lions share are less than an acre. That’s their desired condition, after they treat it. Before treatment, they may be 100% full because they haven’t thinned it yet.
Table 3: Horizontal forest structure (landscape pattern) metrics.
| Landscape structure metric | Set | FRAGSTATS acronym |
|---|---|---|
| Mean patch size | Class | AREA_MN |
| Edge contrast | NA | TECI |
| Patch aggregation | Class | AI |
| Nearest neighbor distance | NA | ENN_MN |
| Patch shape complexity | Class | SHAPE_MN |
| Patch shape complexity | Class | FRAC_MN |
| Patch shape complexity | Patch | FRAC_AM |
| Patch shape complexity | Class | FRAC_CV |
| Patch dispersion | NA | ENN_CV |
| Large patch dominance | Class | LPI |
| Large patch dominance | Patch | AREA_AM |
| Large patch dominance | Patch | CORE_AM |
| Large patch dominance | NA | DCORE_AM |
| Shape and correlation of large patches | NA | GYRATE_AM |
| Shape and correlation of large patches | Patch | SHAPE_AM |
Let’s say we’re interested in generating landscape pattern metrics in a 1-acre moving window for the project area, to develop the data necessary to calibrate the model we:
Figure 3. An example of a simulated forested area with and without (left and center panel, respectively) at a 1-acre scale of analysis. As a result of the presence of shadows (among other factors), estimates of LPIs can be biased high or low, relative to the true value (right panel).
Data from the simulations allow us to fit a model to predict the ‘true’ value of each LPI, i.e., the value we would have obtained if shadows were not present. We modeled the data as
yij = β0 + β1x1, ij + β2x2, i + β3x3, i, ..., βmxm, i + ϵij
where yij is the true value of metric i for sublandscape j. Predictors include the observed value of the metric x1, ij, the proportion of canopy in the sublandscape, x2, j, and the ratio of the proportion of shadow to the proportion of canopy in the sublandscape, x2, j. Parameters β4 − m are tied to the two way interactions between x1 − 3 We then leveraged the relationship between true and observed values to correct each metric for shadow effects, and generated an estimate of our degree of confidence in the reported LPI as the variance of observations around the 1:1 line (?). As expected, shadows have a stronger influence on some variables than others (e.g., AI vs. FRAC_CV).
Figure 4: Scatterplots of the observed vs. true values of each LPI. Perfect estimates would fall along the 1:1 line. All others are either biased high or low.
Figure 5: Values of each LPI against true values after calibration [TODO: report fit statistics.].
[TODOS: apply the model to one of the LPIs and show a before/after. Define functions for metrics we’re currently missing and generate those too. Restart the CCC (V2?) without NAIP.]
Guidance for the acceptable use and appropriate scales of use for each product…. We recommend visual inspection of the classification (and any derived quantities) for class confusion in green meadow areas, or where Pondersoa gives way to scrub brush communities near the rim. Steep terrain…. Changes in canopy cover over time….
All of the code used in developing this analysis is available here.
L. J. Zachmann and B. G. Dickson received support from the the USDA Forest Service, Coconino National Forest (challenge cost share supplemental project agreement 15-CS-11030420-033). Christopher Ray, Valerie Horncastle, and Michael Peters contributed to training (and other ancillary) data development.
Anselin L. 1995. Local indicators of spatial association—LISA. Geographical Analysis 27:93–115. Available from http://onlinelibrary.wiley.com/doi/10.1111/j.1538-4632.1995.tb00338.x/abstract (accessed October 25, 2016).
Breiman L. 2001. Random forests. Machine Learning 45:5–32. Available from http://link.springer.com/10.1023/A:1010933404324 (accessed October 25, 2016).
Conners RW, Trivedi MM, Harlow CA. 1984. Segmentation of a high-resolution urban scene using texture operators. Computer Vision, Graphics, and Image Processing 25:273–310. Available from http://www.sciencedirect.com/science/article/pii/0734189X8490197X (accessed October 25, 2016).
Cushman SA, McGarigal K, Neel MC. 2008. Parsimony in landscape metrics: Strength, universality, and consistency. Ecological Indicators 8:691–703. Available from http://www.sciencedirect.com/science/article/pii/S1470160X07001306 (accessed September 21, 2016).
Farr TG et al. 2007. The shuttle radar topography mission. Reviews of Geophysics 45:RG2004. Available from http://onlinelibrary.wiley.com/doi/10.1029/2005RG000183/abstract (accessed October 25, 2016).
Haralick RM, Shanmugam K, Dinstein I. 1973. Textural features for image classification. IEEE Transactions on Systems, Man, and Cybernetics SMC-3:610–621.
Theobald DM, Stevens DL, White D, Urquhart NS, Olsen AR, Norman JB. 2007. Using GIS to generate spatially balanced random survey designs for natural resource applications. Environmental Management 40:134–146. Available from http://link.springer.com/article/10.1007/s00267-005-0199-x (accessed September 21, 2016).